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NONPARAMETRIC ESTIMATION OF THE PROBABILITY OF 
A LONG DELAY IN THE M/G/1 QUEUE 

BY 

D. P. GAVER 
P. A. JACOBS 

OPERATIONS RESEARCH DEPARTMENT 
• NAVAL POSTGRADUATE SCHOOL 
MONTEREY, CA 93943-5000 

1 . Introduction 

The application of probability theory to a wide variety of congestion 
problems has been well catalogued in many papers. Most of the elegant 
solutions obtained are in somewhat implicit form, being presented as 
functional equations, or, frequently, as integral (Laplace) transforms, 
generating functions, and sometimes as combinations of the above. The 
results obtained naturally appear in terms of component distribution 
functions and stochastic processes (renewal, Poisson etc). Only rarely are 
issues addressed that arise when actual data is to be used as a basis for 
inference from the models; however, see Cox [1965]. 

In this paper we consider the nonparametric estimation of the 
probability of a long customer delay in an M/G/1 system, given a known 
Poisson arrival rate X and observations of independent service times from 
the service distribution, presumed unknown. Although the approach and 
results are given concretely for the M/G/1 system, they apply more widely. 

To be specific, consider a single server system approached by a 
stationary Poisson (X) traffic with X known . Service times, X are 
independent identically distributed with unknown distribution function F^; 
assume XE[X] = p < 1 . Let observations of the service times be all that is 



known about F; denote these by x^. The objective is to supply 

nonparametric estimates, and error assessments thereof, of the probability 
of a long delay experienced by an arriving customer. 

It is well known that if W(t) is the virtual waiting time in the M/G/1 
queue and p < 1 , then the moment generating function 



where A(s) is the moment generating function of a distribution H. If A(s) 
exists for s < s^, s^ > 0, then there will be a smallest real zero s = k > 0 
of the denominator of (1.1) which can be used to show that 



„ r sW 1 , . _ r sW ( t ) 1 

E[e J = lim E[e J 



( 1 . 1 ) 







= (1-p)[1 - pA(s)]“-’ 



P{W > w} - D(k) e 



W -*• 00. 



( 1 . 2 ) 



We will always assume this is the case. 



One way of establishing (1.2) is to introduce 



E[e 



sW 



1 ] 



00 



3 




(1.3) 



0 



into (1.1) and to rewrite in the form 




+ pA(s)H'(s) 



(1.4) 
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which is equivalent to the terminating renewal equation, see Feller [1966] 
p. 362, 



__ w 

F^(w) = P{W > w} = pH(w) P I > w-x} H(dx) (1.5) 

0 



where 

H(w) = I dy. (1.6) 

0 

— 

Following Feller, introduce F^ (w) = e , 0 real and positive, into 

(1.5) to obtain 



^W^^(w) = P H^*^(w) + 



I F^'^(w-x) pe^^H(dx). 
0 



(1.7) 



Choose 0 = < so that 



I pe^^H(dx) = I H^‘^(dx) = 1 
0 0 



( 1 . 8 ) 



yielding a standard renewal equation for F^. From the key renewal theorem, 
it follows that 



lim F (w) 
w 

W^CO 



, . - , s <W 

lim F,,(w)e 
w 

W->oo 



c(k) 

m(<) 



(1.9) 
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where 



and 



c(<) “ P I H(x)dx 

0 



m(ic) = p I xe H(dx). 
0 



Summarizing 



( 1 . 10 ) 



( 1 . 11 ) 



P{W > t} ~ ^ 7 ^ e t 

m(ic) 



( 1 . 12 ) 



where k is the positive solution to the equation 



xe ^[())(e) - 1] = 1 



(1.13) 



with 



<t>(e) = E[e®^]. 



(1.14) 



In section 2, a nonparametric estimate, kc, of < is studied which solves 
equation (1.13) with the empirical moment generating function replacing (p. 
This estimate is related to an estimate studied by Stigler [1971] in the 
context of estimating the probability of extinction for branching processes. 
K is shown to be a consistent estimate of k having a distribution in the 
domain of attraction of a stable law. Under certain conditions, a central 
limit theorem for < can be obtained. 
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In section 3 a nonparametric estimate of 



p(t) 



c( k) 
(k) 



(1.15) 



is given. 

In Section some results of simulation studies of the estimates of < 
and p(t) are presented. 

2. Nonparametric Estimation of the Exponent < of the Probability of a 
Long Wait 

Assume for the remainder of the paper that A=1 is known. Let 
, ..., be the observations of the service times. A non- parametric 
estimate of the moment generating function of the service time distribution 
is 



.. . n 0x . 

<()(e) = - I e \ ( 2 . 1 ) 

" i = 1 



The sample equivalent to equation (1.13) is thus 



1 = e ^[i(e)-1j. (2.2) 



At 6 = 0, the RHS of (2.2) is x which is less than 1 if 

X +X +...+X 

X = ^ < 1; the data can only be analyzed for a stationary 



model under this assumption, which will be made in what follows. Further, 
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(2,2) has a unique solution, k, which is an estimate of k . A three-term 
Taylor's expansion of the RHS of (2.2) about 0=0 gives an initial estimate 

= 2[1 - x] (7^)"’ (2.3) 

”k 1 k k th 

where x = - [x, + ..• + x ], the k sample moment. Since 
n 1 n 

s'' [4>(6) ^ 1] ^ X + T? 6 x^, ic„ is an upper bound for k . Equation (2.2) can 

d H 

be solved via search or Newton-Raphson iteration to obtain the estimate k. 

We will now present asymptotic results for the distribution of k as the 
sample size n ■* ®. By assumption, if X is a random variable having the 
service time distribution, then E[e ] < ® for some B' > k > 0. Thus, 
there is B > k such that for all 0 < b S B 

E[x"e^^] < 00. (2.4) 

It follows from the monotonicity of the RHS of (2.2) that 

Pi; > B) . p|i<£^ - I 

D D D 

= p{i(B) - 4.(B) < B[1 (2.5) 

Since <j)(0) is a monotone function and B > k, [1 - B ^[(|)(B) -1]] is negative. 
Thus, by the strong law of large numbers 
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( 2 . 6 ) 



lim P{ K > B} = 0. 

n^co 



Let 

f(6) = 1 - e"’[4>(e) - 1]. 



(2.7) 



Expand f(<) in Taylor's series about the solution < of (1.13)- Since 
f(<) = 0, 

0 = T(k) + (k-k) f'(Kr) + (k-k)^ f''(0K) (2.8) 



for some 9. Thus 



K - K = -f(K)[f'(K) + 1 ( k - k ) f’'(6K)] ^ (2.9) 

E[f(K)] = 0; 

* 1 k'Y k-Y 

E[f'(K)] = E[-<Xe + e - 1] = 6 < 0. (2.10) 



By the strong law of large numbers 



lim }{ k ) = 0 (2.11) 

n-^oo 

and 

lim f'(ic) = 6 < 0 (2.12) 

n->-“ 
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with probability 1. It follows from (2.6), (2.9), (2.11) and (2.12) that k 
converges in probability to k as n Thus, k is a consistent estimate of 

K. 



?krY 

If E[e ] < then 

Var[f(<)] = - -4 Var[e^^] = - 
n n 

By the central limit theorem 



(2.13) 



C-f(K)] Cf'(K) + ^ (k-k) f’*(6K)] ’ (2.U) 



is asymptotically standard normal as n 

2<X kX 

If E[e ] = “, then the distribution G of e is in the domain of 

attraction of a stable law with index 1 < a < 2, see Feller [1966]. Let 

{a^} be a sequence of numbers such that 




G(dx) ^ 1 



as n ^ 



(2.15) 



The normalized random variable 



A . n kX , .A y. . 

^ [k - k] = I e h[f'(<) + ^(k - k) f'(0K)] 

a an., 2 

n n 1=1 



( 2 . 16 ) 
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is in the domain of a stable law with index a, where 1 < a < 2 is such that 



I 



0 



x^G(dx) ~ ^ 

To summar i ze , 



L(y) where L is a slowly varying function, 
we have the following result. 



PROPOSITION a. < is a consistent estimate of k. 

2kX 

b. If E[e 3 < ”, then the distribution of <-< is 

asymptotically normal. 

2<X 

c. If E[e ] = ”, then kt-k is in the domain of 
attraction of a stable law. 

We will suppose for the remainder of this section that the service time 
distribution is exponential with mean ~ In this case, (since A=1 by 

assumption) , 



< = u - 1 : (2.17) 

i^Y P - 1 

Var[e ] = y(y-l) (2-y) ; (2.18) 

and 

ECf'(K)] = -1. (2.19) 

From the central limit theorem (2.1^1) 

nVar(<) - = y(2-y)"’ = (2.20) 
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2 1/2 -1 

Therefore, if g(x) = (1-x ) - 1 - sin (-x), then 



nVar[g(K)] ~ [g'(i<)] nVar[i<] ~ 1. 



Thus, g is an asymptotically variance-stabilizing transform of k. The 
Simpler related transf ormation ln(1+<) is used in the simulations of < 
reported in Section 4. * 

3. A Nonparametric Estimate of P{W > t} 

The asymptotic analytic result is 




(3.1) 



as t ® 



where 



00 



00 




(3.2) 



0 y 
00 

= J 1^ P{Xedx} + ^ 



00 



f r 

2 J Ce 



00 



KX 




0 



0 



and 



00 




(3.3) 



0 
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= ^ (j Kxe"" P{Xedx} - j - 1) 

0 0 

CD 00 

= ^ I x^ P{Xedx} * ^ ^ j P{Xedx} 

0 0 
00 

+ I (kx - 1) R(x) P{Xedx} 



where 



r>/ \ r <x ; I 

R(x) = [e - 1 - <x 5 — J 



(<x)' 



An estimate of P{W > t} is 



^ . s c ( k) -Kt 
p(t) = -s— s — e 



m(<) 



where k is the positive solution to equation (2.2) 



c(ic) = + T— R^ ( k) ; 

m(<) = + ^ < x’ + R^Ck); 

Rl(<) = - J (e " - 1 



.S (kx. ) 

- KXj - -ji— ), 



P{ Xedx} I 



( 3 . 4 ) 



( 3 . 5 ) 



( 3 . 6 ) 



( 3 . 7 ) 



( 3 . 8 ) 



-n- 



1 



(3.9) 



H^U) = (- R^(k)) + < ^ I X. 

i = 1 



KX . 
1 



- <x. 
1 



(kx ) 

). 



The forms of c(k) and m(K) are chosen for the numerical stability of the 
ratio c(k) m(K) ^ . 

In the remainder of this section we will assume the service time 
distribution is exponential with mean “ > ^ x < 1. In this M/M/1 queue 
case, it is well known that 



p(t) = exp{-(u“1)t} = (1 + k) ^ e (3*10) 

f r om ( 2 . 1 6 ) . 

We will now motivate a transformation of p(t) which is used in the 

Y 

simulation studies. Let Y = ln(1+K), then k = e - 1 and 



p(t) = exp{-Y - [e^ - l]t}. 



(3.11) 



Let 



h(Y) = In p(t) = - Y - [e - 1]t 



and Y = ln(1 + k) . A Taylor's expansion yields 



nVarCY] = nVar[ln(1 + <)] 

— = [1-(e - 1 ) ] 



1 -ic 



( 3 . 12 ) 



(3.13) 
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from (2.20). Hence, 



nVar[h(Y)] = [-1 -e\]^ [1 - (e"^ - 1)^] ’ 

= C- 1 - t - (e^ - 1)t]^[1 - (e"^ - 1)^]“''. 



( 3 . 1^0 



It follows from the definition of h(Y) that 



nVarCh(Y)] = [-(H-t) + h(Y) + Y]^{1 - (e^ - 1)^) ’ 



(3.15) 




nVar[h(Y)] = [-h(Y) + t]^ 



(3.16) 



which suggests that the transformation ln[t - In p(t)] will tend to 
stabilize the variance of p(t), at least for exponential service time 
distributions . 

Simulation Results 

In this section results are presented of a simulation experiment to 
study small sample behavior of the estimates of < and P{W > t}. All 
simulations were carried out on an IBM 3033AP computer at the Naval 
Postgraduate School using the LLRANDOM II random number generating package; 
Lewis and Uribe [1981]. In each replication 50 exponential service times 
are generated. If the mean of the service times is larger than 1, equation 
(2.2) has no positive solution. In this case, another 50 service times are 
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generated. The estimates, < of (2.2) and p(t) of equation (3.5) are 
computed for each replication. 

Standard errors of the estimates are estimated in two ways. One is the 
appropriate asymptotic variance expression. The other is the Jackknife 
procedure. 

The jackknife is a procedure originally introduced by Quenouille for 

bias reduction [1956], and adapted by Tukey [1958] to obtain approximate 

confidence intervals. Suppose interest is on a parameter 0 (e.g. k or p(t)) 

that is estimated by 0 using a complex calculation from data x^,...,x^. The 

idea is that of assessing variability by recomputing 0 after removing 

independent subgroups of data of equal size, and then using the recomputed 0 

values to estimate a variance which is in turn applied to state a standard 

error or a two-sided confidence interval that contains the true 0 with 

specified confidence. A few more details follow; for more, see Efron [1982] 

and his more recent work, or Hosteller and Tukey [1977]. The actual 

calculation involves splitting the n data points into g disjoint groups of 

size m; n=mg. In our simulations g is always 10. Then calculate 

j=1,2,...,g: the estimate of 0 based on a reduced data set that omits the 
t h 

j group. In the simulations, the first group is the first five 
(unordered) service times; the second group is the second five service 
times, etc. 

Now Tukey computes pseudo-values 



Yj = ge 



(g-1) 0 



(-J) 



-in- 



which are treated as independent. Tukey recommends referring the mean of 
the pseudo-values y to Student's t with g-1 degrees of freedom to obtain 
confidence limits. 

In the jackknife procedure to estimate k, it is sometimes the case that 
a positive solution of equation (2.2) does not exist for the data set that 
omits the j^*^ subgroup because the mean of this reduced data set is larger 
than 1. In this case, is set equal to the smallest of the k^^'s that 

could be computed for the sample considered.. 

Similiarly, in a jackknife procedure to estimate p(t), it is sometimes 

^ /S 

the case that either does not exist or p^_j^(t) exceeds 1 for the 

reduced data set that omits the subgroup. In this case is set 

equal to the largest of those that could be computed and were less 

than 1 , 



4. 1 Results of a Simulation Experiment to Estimate < 

In this subsection results are given of a simulation experiment to 

estimate kc. For each replication, three different estimates of ic were 

computed. Estimate I, , is computed by numerically solving (2.2) by 

search starting with the initial value of (2,3)* Estimate II, ic 

H 11 

is obtained by jackknifing k using ten subgroups; the mean of the pseudo- 
values is the estimate. Estimate III is obtained by jackknifing ln(l+K) 
using ten subgroups; the inverse transform of the mean of the pseudo- values 

- 1 is used as the estimate of c. In the cases of Estimates II and III, 
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if the estimate is negative, it is set equal to 0. For each estimate, the 



average bias 



B 



, 500 

— y 

i=i 






and the relative mean square error 



Rel MSE 



1 

500 



500 <. -< 2 
I [-^] 

i=1 ^ 



are computed where 

ic = p-1 

in the case of exponential service times with mean ^ 1 • 

Results of the simulation appear in Table 1. All of the estimates have 
about the same relative mean square error. Jackknife estimates II and III 

have smaller bias than the straightforward estimate I. Jackknifing < itself 

^ 1 
rather than ln(<+1) gives the smallest bias. As — increases all the 

y 

estimates have increased relative mean square error. 

A simulation study was conducted to compare the performance of 
different confidence interval procedures. For each replication three 80^ 
confidence intervals were constructed. Interval procedure I is a normal 
confidence interval which uses the straightforward estimate of k, k:^, as the 
point estimate; the estimate of the variance is the data version of the 
asymptotic variance in the central limit theorem (2.14); 
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i 



2 



(^. 1 ) 



where 



^ I e ^ 



i 



(^. 2 ) 



<x. e 
1 



(^. 3 ) 



i 



with K = Kj in this case. The 80^-point of the normal distribution is used 
to construct the interval. Limits that are negative are set equal to 0. 

Confidence interval procedure II is a jackknife confidence interval 
which Jackknifes < and uses the 80% point of the student t-distribution with 
9 degrees of freedom. Limits that are negative are set equal to 0. 

Confidence interval procedure III is a jackknife confidence interval 
which jackknifes ln(<+1) and uses the 80% point of the student 
t-distribution with 9 degrees of freedom to give a confidence interval for 
ln(<+1). The inverse transformation of the endpoints of the interval gives 
an interval for <; limits that are negative are set equal to 0. 

Results of the confidence interval simulation appear in Table 2. 
Reported are the number of 500 intervals that cover the true value of 
K, (C); the number of the 500 intervals such that the entire interval lies 
below K, (B); and the number of the 500 intervals such that the entire 
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interval lies above the true value, (H). The average length of the 
confidence intervals is also given. 

The number of intervals that cover the true value of k for procedure 

III is within .80 ± (1.96)X (.2) (.8) = [.765, .835]. All but one case 

duu 

of the normal confidence intervals of procedure I are outside this range. 
All but one case using confidence interval procedure II are inside this 
range. The average width of confidence intervals for procedures II and III 
are about the same. Thus, although the jackknife estimate is a little 

more biased than K^j,the coverage of the Jackknife confidence interval for 
appears to be somewhat better than for Jackknife confidence interval 
procedure II. 
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TABLE 1 

Bias and Mean Relative Square Error 
for Estimates of k 



Distribution 

Estimate 


Exponential 

- = 0.6 
M 

< = .6667 


Exponential 

- = 0.7 
u 

< = .4286 


Exponential 

- = 0.8 
u 

K = .2500 


Exponential 

- = 0.9 
u 

K = .1111 




B Rel MSE 


B Rel MSE 


B Rel MSE 


B Rel MSE 




(S.E) 


(S.E. ) 


(S.E) 


(S.E. ) 


(S.E) 


(S.E.) 


(S.E) (S.E.) 




(-) 




(f) 




(f) 




(f) 






K 
















I 


,1147 


.2350 


.0557 


.3292 


.0057 


.6026 


.0865 


2.400 




(.014) 


(.021 ) 


(.011) 


(.029) 


( .008) 


(.054) 


( .0067) 


(.2134) 




[.172] 




[.130] 




[.023] 




[.779] 




II 


.0219 


.2256 


.0067 


.3095 


.01 02 


.5382 


.0469 


1.785 




( .014) 


( .019) 


(.011) 


( .025) 


( .008) 


(.046) 


( .0063) 


(.161 ) 




[.033] 




[.016] 




[.041 ] 




[.422] 




III 


.0533 


.221 2 


.0160 


.3073 


.0268 


.5576 


.0608 


1.997 




(.014) 


( .019) 


(.01 1) 


( .027) 


( .008) 


(.049) 


( .008) 


(.177) 




[.080] 




[.037] 




[.1072] 




[.5472] 
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Table 2 



Confidence Intervals for < 



Procedure 

Distribution 


Coverage 


Width 


B 

% 


I 

(C) 

(%) 


. H 
% 


B 

% 


II 

(C) 

(%) 


H 

% 


B 

% 


III 

(C) 

i%) 


H 

% 


I 

AV 

(S.E.) 


II 11 

AV AV 

1 (S.E.) (S. 


Exponential 


























- = 0.6 
|i 


45 


( 363 ) 


92 


38 


( 408 ) 


54 


29 


( 395 ) 


76 


.6402 


.8089 


.77 


K = .6667 


9 


( 72 . 6 ) 


18.4 


7.6 


(81 . 6 ) 


10.8 


5.8 


( 79 ) 


15.2 


(. 005 ) 


(. 0124 ) 


(.01 


Exponential 


























- = 0.7 

P 


49 


( 380 ) 


71 


50 


( 409 ) 


41 


42 


( 399 ) 


59 


.5279 


.6255 


.61 


K = .4286 


9.8 


( 76 ) 


14.2 


10 


( 81 . 8 ) 


8.2 


8.4 


( 79 . 8 ) 


11.8 


(. 004 ) 


(. 009 ) 


(.00 


Exponential 


























- = 0.8 
]i 


28 


( 407 ) 


65 


38 


( 420 ) 


42 


33 


( 413 ) 


54 


.4441 


.4831 


.45 


K = .2500 


5.6 


( 81 . 4 ) 


13 


7.6 


( 84 ) 


8.4 


6.6 


( 82 . 6 ) 


10.8 


(. 005 ) 


( . 008 ) 


(.00 


Exponential 


























- - 0.9 

p 


1 


( 429 ) 


70 


54 


( 408 ) 


38 


52 


( 392 ) 


56 


.3733 


.3763 


.39 


K = .1111 


. 002 ( 85 . 8 ) 


14 


10.8 


( 81 . 6 ) 


7.6 


10.4 


( 78 . 4 ) 11 .2 


(. 005 ) 


(. 009 ) 


(.0 



M . 2 Simulation Results for Estimating p(t) 

In this subsection, results are given of a simulation experiment to 
estimate p(t). For each replication three estimates of p(t) are computed. 
Estimate I, Pj(t) is computed for each replication using formulas (2.2), 
(3-5) - (3.9). If Pj(t) exceeds 1, it is set equal to 1. 
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There are (at least two) possible ways to implement a jackknife 
procedure to estimate p(t). Estimate II is obtained by jackknifing ln(<+l); 
an estimate of k is obtained by the inverse transformation of the mean of 
the pseudo-values y^^^, 




1 ; 



Kjj is used in formulas (3.5) - (3.9) to obtain the estimate p^j(t). If 
Pjj(t) exceeds 1, it is set equal to 1. 

Estimate III is obtained by jackknifing ln[t - In p(t)]; if p(t) 
exceeds 1 for a reduced data set that omits the j^*^ subgroup, the estimate 
of ln[t - In p(t)] for that reduced data set is put equal to the smallest 
estimate that could be computed from the other reduced data sets. An 
estimate of p(t) is obtained by inverse transformation of the average of the 
pseudo-values y, no: 

Lr K 



t "^LPR 
Plll(t) = e e 



If Pjjj(t) >1, it is replaced by 1. 



For each estimate, the average bias 

. 500 n 

® “ 500 ^ tPi(t) - p(t)l 

i = 1 



(4.iJ) 



and the relative mean square error 



2 



Rel MSE - ^ 



500 



(t ) - p(t ) 
p(t) 



(^. 5 ) 



are computed where 



p(t) 




1)t 



(M.6) 



The results appear in Table 3* 



TABLE 3 

Mean Bias and Mean Relative Square Errors 
for Estimates of p(t). 



Distribution 


Exponential 

- =.6, T=1 
U 

p(1) = .3081 


Exponential 
1 -.7, T.2 
p(2) = .2971 


Exponential 
1 -.8, T-3 

p(3) = .3779 


Estimate 


B Rel MSE 


B Rel MSE 


B Rel MSE 




(S.E) (S.E.) 


(S.E) (S.E.) 


(S.E) (S.E.) 


I 


.00^1 .138 


.024 .328 


.005 .326 




(.005) (.009) 


(.008) (.026) 


(.010) (.020) 


II 


.044 .216 


.060 .409 


.054 .436 




(.006) (.020) 


(.008) (.032) 


(.011) (.031) 


III 


.0138 .137 


.040 .354 


.034 .397 




(.005) (.009) 


(.008) (.030) 


(.011) (.029) 
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The entries for estimate II in Table 3 suggest that jackknifing < and 
then computing the estimate of the probability using the jackknifed estimate 



of K increases the bias and relative mean square error of the estimate over 



the jackknifing InCt - In p(t)] gives about the same relative mean square 
error as the original estimate Pj(t) but increases the bias somewhat. 

A simulation experiment was conducted to investigate the performance of 
confidence intervals obtained by jackknifing ln[t - In p(t)]. For each 
simulation replication, the average and variance of the pseudo-values are 
computed; 



where 1.383 is the 80% point for a Student's t-distribution with 9 degrees 
of freedom. The limits of the interval are inversely transformed to give a 
confidence interval for p(t). 



Pj(t) the straight-f orward estimate. The entries for estimate III suggest 





An Q0% confidence interval for ln[t - In p(t)] is 
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PL = exp{t - 



and 

PV = exp {t - . 

If a confidence limit exceeds 1, the limit is set equal to 1. 

Results of the confidence interval simulation appear in Table 4. 
Reported are the number of the 500 intervals that cover the true value 
p(t) (C); the number of the 500 intervals for which the upper limit PU is 
below p(t) (B); and the number of the 500 intervals such that the lower 
limit PL is above the true value (H). The average width of the confidence 
interval is also given. 
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Table 4. 

80? Confidence Intervals for p(t) 







Coverage 




Average Width 




B 

(?) 


C 

(?) 


H 

(?) 


(standard error) 


Exponential 










- =.6, T=1 

p 


56 


41 5 


29 


.325 


p(1) = .3081 


(11.2) 


(83) 


(5.8) 


( .006) 


Exponential 










- =.7, T=2 


52 


408 


40 


.477 


U 

p(2) = .2971 


(10.4) 


(81 .6) 


(8) 


(.009) 


Exponential 

1 -8, T-3 


54 


412 


34 


.597 


p(3) = .3779 


(10.8) 


(82.4) 


(6.8) 


(.010) 



The coverage of the confidence intervals is within 

.80 ± (1.96) ‘^ 500 ^^ = *80 ± *035 = [.765, .835]. The width of the 

interval increases as — increases. 

P 

4.3 Simulation Results for Estimating p(t) for Mixed Exponential Service 
Times 

In this subsection, results are given of a simulation experiment to 
estimate p(t) = P{W>t} in the case of mixed exponential service times. For 
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each replication 50 random numbers are generated from a mixed exponential 
distribution with 



1 1 

P{S>t} = ^ ^ , t>0. 



The estimate ln[t - Inp(t)] is jackknifed with ten subgroups as before where 
p(t) is given by (3.5) and 80? confidence intervals are constructed. Each 
simulation has 500 replications. 

The asymptotic distribution of the virtual waiting time, W, can be 
found by inverting the moment generating function (1.1); it is a mixture of 
two exponential random variables. 

Table 5 reports results of three simulations. 



Case I : 



1_ 

E(S) 



.30 



= 0 . 6 , 



- = . 90 , T = 1 , 

P{W>T} = .3^011; 



Case II: 



1_ 

^*1 

E(S) 



.35 



0.7, 



- = 1.05, T = 2, 

^*2 

P{W>T} = .3^59; 



Case III; 



i = .l;0 - = 1 .20, T = 3, 

Ul U2 

E(S) = 0.8, P{W>T} = .ii333; 
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For each simulation, the mean bias and relative mean square error 
and (4.5) are computed with p(t) = P{W>t}. Confidence interval coverage and 
average widths are also computed. 

Both cases II and III each had one replication for which using all the 
data to compute the estimate p(t) of (3.5) resulted in a value larger than 
1 ; these two replications are not counted in the summary statistics in Table 

5. 

Comparison of the mean biases and mean relative square errors of Table 
5 with those of estimate III in Table 3 shows that they are about the same 
for both service time distributions. The coverage of the confidence 
intervals in Table 5 is again within [.765, .835]. However, the average 
length of the confidence intervals for the mixed exponential cases are 
larger than those reported in Table 4 for the exponential distributions. 
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Table 5 



Results for Estimates of P{W>t} 
for an M/G/1 Queue with Mixed 
Exponential Service Times 



Distribution 

Case 


Estimate 
Bias Rel MSE 
(S.E) (S.E.) 


Conf idence 
Coverage 
B (C) H 

% (?) ? 


Intervals 

Width 

Average 

(S.E.) 


I 


.024 


.172 


58 


(407) 


35 


.392 




( .006) 


(.013) 


11.6 


(81.4) 


7 


(.007) 


II 


.040 


• .384 


63 


(408) 


28 


.500 




(.009) 


(.028) 


12.6 


(81 .4) 


5.6 


( .010) 


III 


.022 


.350 


54 


(411) 


34 


.637 




(.011) 


(.021 ) 


10.8 


(82.4) 


6.8 


(.011) 



4.i| Simulation Results for Estimating P{W>t} for Gamma Distributed Service 

Times 

In this subsection, results are given of a simulation experiment to 
•estimate p(t) = P{W>t} in the case of gamma service times. Each experiment 
has 500 replications. For each replication, 50 service times are generated 
having the distribution of the sum of two exponential random variables each 
having mean The estimate ln[t - In p(t)] is jackknifed with ten 
subgroups as before where p(t) is given by (3.5) and 80? confidence 
intervals are constructed. 
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The asymptotic distribution of the virtual waiting time, W, can be 
found by inverting the moment generating function (1.1); it is a mixture of 
two exponential random variables. 

Table 6 reports results of three simulations. 



Case I : 




E(S) = 0.6, P{W>T} = .2508; 



Case II : 




T = 2, 



E(S) = 0.7, P{W>T} = .2232; 



Case III: 



- = .40 T = 3, 

M 

E(S) = 0.8, P{W>T} = .2950; 
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Table 6 



Results for Estimates of P{W>t} 
for an M/G/1 Queue with Gamma 
Service Times 



Distribution 

Case 


Estimate 
Bias Rel MSE 
(S.E) (S.E.) 


Confidence 
Coverage 
B (C) H 

% (%) % 


Intervals 

Width 

Average 

(S.E.) 


I 


..007 


.110 


61 


(ilOl) 


38 


.227 




(.00^0 1 


(.007) 


1 2.2 


(80.2) 


7.6 


(.003) 


II 


.0167 


.299 


61 


(397) 


H2 


OO 

C\J 

on 

• 




(.005) 1 


(.026) 


12.2 


(79.^) 


8.1) 


( .007) 


III 


.0251 


.390 


62 


(i107) 


31 


.1190 




(.008) 1 


(.03^) 


1 2.H 


(81. H) 


6.2 


(.010) 



The mean bias and relative mean square error are about the same as for 
the exponential and mixed exponential service time distributions. The 
confidence interval coverage is once again within [.765. .835]. The average 
width of the confidence intervals is smaller than the widths in the 
exponential and mixed exponential cases. This is to be expected since the 
gamma distribution has a shorter tail than the exponential or mixed 
exponential distribution. 
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